#define EPSILON 0.0001f
#define MAX_STEP 1000

uniform sampler2D u_TransferFunctionTex;
uniform sampler3D u_VolumeTex;

uniform bool u_UseGradient;
uniform float u_SamplingFactor;

ivec3 Dimension = textureSize(u_VolumeTex, 0);
vec3 DimensionInterval = 1.0f / vec3(Dimension);

//TODO: precompute to accelerate?
// Gradient on-the-fly using tetrahedron technique http://iquilezles.org/www/articles/normalsSDF/normalsSDF.htm
float getGradient(vec3 vPosition, vec3 vDimensionInterval)
{
	float GradientMagnitudeModifier = 1.0f;
	ivec2 k = ivec2(1, -1);
	vec3 GradientDirection = (k.xyy * texture(u_VolumeTex, vPosition + vDimensionInterval * k.xyy).x +
								k.yyx * texture(u_VolumeTex, vPosition + vDimensionInterval * k.yyx).x +
								k.yxy * texture(u_VolumeTex, vPosition + vDimensionInterval * k.yxy).x +
								k.xxx * texture(u_VolumeTex, vPosition + vDimensionInterval * k.xxx).x) * 0.25f;

	return clamp(length(GradientDirection) * GradientMagnitudeModifier, 0.0f, 1.0f);
}

bool equals(float x, float y) {
	return abs(x - y) < EPSILON;
}

vec3 modifyVecNE2Zero(vec3 Vec) {
	if (equals(Vec.x, 0.0f)) Vec.x += EPSILON;
	if (equals(Vec.y, 0.0f)) Vec.y += EPSILON;
	if (equals(Vec.z, 0.0f)) Vec.z += EPSILON;
												  
	return Vec;
}

vec3 computeStepSize(vec3 vRayEntry, vec3 vRayExit, vec3 vDirection, out int vNumSteps = 0) 
{
	float RayVolumeDistance = distance(vRayEntry, vRayExit);
	
	float MaxDimension = max(Dimension.x, max(Dimension.y, Dimension.z));
	vNumSteps = int(ceil(MaxDimension * RayVolumeDistance * u_SamplingFactor));

	float Denom = float(vNumSteps) - 1.0f;
	if (equals(vNumSteps, 0.0f)) Denom += EPSILON;
	vec3 StepSize = vDirection * RayVolumeDistance / Denom;

	return StepSize;
}

vec4 accumulateColor(vec3 vOrigin, int vStepStart, int vStepEnd, vec3 StepSize) 
{
	vec4 Color = vec4(0.0f);
	float SampleFrequency = 1.0f / u_SamplingFactor;

	float Gradient = 1.0f;
	int NumStep = vStepEnd - vStepStart;
	for (int i = 0; i < min(NumStep, MAX_STEP); ++i)
    {
		vec3 SamplePos = vOrigin + float(i + vStepStart) * StepSize;
		float Intensity = texture(u_VolumeTex, SamplePos).r;
		if(u_UseGradient)
			Gradient = getGradient(SamplePos, DimensionInterval);
		
		vec4 SampleColor = texture(u_TransferFunctionTex, vec2(Intensity, Gradient));
		SampleColor.a = clamp(1.0f - pow(1.0f - SampleColor.a, SampleFrequency), 0.0f, 1.0f);

		// from back to front
		//vec4 SampleColor = texture(u_TransferFunctionTex, vec2(Intensity, Gradient));
		//Color = SampleColor * SampleColor.a + (1 - SampleColor.a) * Color;
		
		// from front to back
		SampleColor.xyz *= SampleColor.a;
		Color = Color + (1.0f - Color.a) * SampleColor;

		if (Color.a >= 1.0f) break;
    }

	return Color;
}